module modsystem
    use modzrodlo
    use ifport
    use modinip
    implicit none
    private

     ! -----------------------------------------------------------------
     ! Zmienne systemowe
     ! -----------------------------------------------------------------
    integer :: no_zrodel                ! liczba zrodel w ukladzie
    integer :: TRANS_MAXN               ! liczba oczek brane do rozwiazania ukladu rownan
    integer :: TRANS_COMPOMENT_MAXN     ! liczba oczek siatki pojedynczego kompenentu
    integer :: MATASIZE
    integer :: Nx,Ny          ! wymiar ukladu
    integer :: no_components  ! liczba skladowych funkcji falowej
    double precision :: DX,Ef,Bz
    doubleprecision  :: TRANS_R ! prawdopodobienstwo odbicia
    doubleprecision  :: TRANS_T ! prawdopodobienstwo przejscia
    doubleprecision  :: TRANS_R_COMPONENTS(500) ! prawdopodobienstwo odbicia dla skladowych f.f.
    doubleprecision  :: TRANS_T_COMPONENTS(500) ! prawdopodobienstwo przejscia dla skladowych f.f.
    double precision :: decoh_amp               ! amplituda dekoherencji

    type(czrodlo),dimension(:),allocatable        :: zrodla ! TABLICA Z OBIEKTAMI ZRODLA
    integer,dimension(:,:), allocatable           :: GFLAGS ! ZAWIERA FLAGI UKLADU (DIRICHLET,WEJSCIA,POZA)
    integer,dimension(:,:), allocatable           :: ZFLAGS ! ZAWIERA FLAGI WEJSC WARTOSC FLAGI OZNACZA NUMER WEJSCIA - 1,2, ...
    integer,dimension(:,:), allocatable           :: GINDEX ! INDEKSUJE GEOMETRIE UKLADU (-1) OZNACZA POZA UKLADEM
    double precision,dimension(:,:), allocatable  :: UTOTAL ! MACIERZ POTENCJALU EFEKTYWNEGO
    complex*16,dimension(:),allocatable           :: CMATA  ! GLOWNA MACIERZ PROGRAMU W FORMACIE (ROW,COL,VALS), TUTAJ TYLKO VALS
    integer,dimension(:,:),allocatable            :: IDXA   ! INDEKSY MACIERZY (ROW,COL)
    complex*16,dimension(:),allocatable           :: VPHI   ! SZUKANA FUNKCJA FALOWA, PRZELICZANA POTEM NA LDOS
    complex*16,dimension(:,:,:),allocatable       ::  PHI   ! FUNKCJA FALOWA ZAPISANA NA DWUWYMIAROWEJ MACIERZY
    complex*16,dimension(:,:,:),allocatable       :: W_kl   ! tablica elementow laczacych skladowe f.f. (i,j,(x,y)) - i,j polaczenie pomiedzy
                                                            ! i-ta, a j-ta skladowa. (x,y) -> indeks wezla siatki
    ! rodzaje zapisu do pliku
    ENUM,BIND(C)
        ENUMERATOR :: ZAPISZ_FLAGI      = 0
        ENUMERATOR :: ZAPISZ_POTENCJAL  = 1
        ENUMERATOR :: ZAPISZ_KONTUR     = 2
        ENUMERATOR :: ZAPISZ_INDEKSY    = 3
        ENUMERATOR :: ZAPISZ_PHI        = 4
        ENUMERATOR :: ZAPISZ_WKL        = 5
    END ENUM

    public :: ZAPISZ_FLAGI , ZAPISZ_POTENCJAL , ZAPISZ_KONTUR , ZAPISZ_INDEKSY , ZAPISZ_PHI , ZAPISZ_WKL
    public :: zrodla,UTOTAL,GFLAGS,GINDEX,PHI
    public :: system_inicjalizacja   , system_zwalnienie_pamieci , system_inicjalizacja_ukladu
    public :: system_zapisz_do_pliku , system_rozwiaz_problem
    public :: system_gauss
    public :: TRANS_T,TRANS_R
    contains


    ! =========================================================
    !               INICJALIZACJA SYSTEMU
    ! =========================================================
    ! PROCEDURA INICJALIZACYJNA WYWOLYWANA JAKO PIERWSZA
    ! PNX,PNY       - ROZMIAR PROBLEMU, WYMIAR MACIERZY GINDEX, GFLAGS
    ! pLiczbaZrodel - LICZBA WSZYTKICH ZRODEL W UKLADZIE
    ! pDX           - Podajemy krok siatki w [nm]
    ! pNo_components- podajemy liczbe skladowych posiadanych przez f.f
    ! ----------------------------------------------------------
    subroutine system_inicjalizacja(pnx,pny,pLiczbaZrodel,pdx,pNo_components)
        integer,intent(in)             :: pnx,pny
        integer,intent(in)             :: pLiczbaZrodel
        double precision, intent(in)   :: pdx
        integer,intent(in)             :: pNo_components
        if(TRANS_DEBUG==.true.) then
            print*,"System: Inicjalizacja:";
            print*,"    nx          =",pnx
            print*,"    ny          =",pny
            print*,"    liczbaZrodel=",pLiczbaZrodel
        endif
        nx            = pnx
        ny            = pny
        dx            = pdx*L2LR
        no_components = pNo_components
        ! alokacja pamieci
        allocate( GFLAGS(nx,ny))
        allocate( ZFLAGS(nx,ny))
        allocate( GINDEX(nx,ny))
        allocate( UTOTAL(nx,ny))
        allocate( PHI(nx,ny,no_components))
        allocate(W_kl(no_components,no_components,nx*ny))
        GFLAGS = B_NORMAL
        ZFLAGS = 0
        GINDEX = 0
        UTOTAL = 0
        ! Tworzymy uklad...
        no_zrodel = pLiczbaZrodel
        allocate(zrodla(no_zrodel))
        ! okreslanie rzeczywistego rozmiaru problemu, punkty o indeksach -1 leza poza ukladem, nie beda
        ! potrzebne wiec zatem do obliczen.
        TRANS_MAXN = 0;

        ! -------------------------------------------------- !
        call getDoubleValue("Dane","decoh_amp",decoh_amp);
        ! -------------------------------------------------- !

    end subroutine system_inicjalizacja


    ! =========================================================
    !
    !           CZYSZCZENIE ZAALOKOWANEJ PAMIECI
    !       funkcja wywolywana jest na samym koncu po wykonaniu
    !       obliczen
    ! =========================================================

    subroutine system_zwalnienie_pamieci()
        integer :: i
        print*,"System: Zwalnianie pamieci"
        if(allocated(GFLAGS))   deallocate(GFLAGS)
        if(allocated(ZFLAGS))   deallocate(ZFLAGS)
        if(allocated(GINDEX))   deallocate(GINDEX)
        if(allocated(UTOTAL))   deallocate(UTOTAL)
        if(allocated(PHI))      deallocate(PHI)
        if(allocated(W_kl))     deallocate(W_kl)
        do i = 1 , no_zrodel
            call zrodla(i)%zrodlo_zwolnij_pamiec()
        enddo
        if(allocated(zrodla)) deallocate(zrodla)
    end subroutine system_zwalnienie_pamieci

    ! ------------------------------------------------------------
    ! Funkcja tworzy tablice flag oraz gindex na podstawie
    ! zmodyfikowanej w mainie tablicy gflags
    ! Podajemy:
    ! in_len - dlugosc na jaka od wejsc zostanie utworzona flaga B_NORMAL
    ! smooth_rad - promien wygladzania ostrych katow w ukladzie
    ! smooth_iter - liczba iteracji wygladzania
    ! ------------------------------------------------------------
    subroutine system_inicjalizacja_ukladu(in_len,smooth_rad,smooth_iter)
        integer ,intent(in) :: in_len,smooth_rad,smooth_iter
        integer :: soft_iter , iii , jjj , norm_iter , empty_iter
        integer :: nrz , i , j , ni, nj , mi , mj , iter

        print * , "System: Tworzenie tablic GINDEX oraz GFLAGS"
        ! dodawanie zrodej wejsciowych do warunkow brzegowych
        do nrz = 1 , no_zrodel
            do i = 2 , zrodla(nrz)%N - 1
                ni = zrodla(nrz)%polozenia(i,1)
                nj = zrodla(nrz)%polozenia(i,2)
                GFLAGS(ni,nj) = B_WEJSCIE
                ZFLAGS(ni,nj) = nrz ! przypisujemy fladze numer wejscia
            enddo
            if(zrodla(nrz)%bKierunek)then
                ni = zrodla(nrz)%polozenia(1,1) + 1
                nj = zrodla(nrz)%polozenia(1,2)
                mi = ni + in_len - 1
                mj = zrodla(nrz)%polozenia(zrodla(nrz)%N,2)
                if(nj==1) nj = 2
                if(mj==NY)mj = NY-1
                GFLAGS(ni:mi,nj:mj) = B_NORMAL
            else
                ni = zrodla(nrz)%polozenia(1,1) - 1
                nj = zrodla(nrz)%polozenia(1,2)
                mi = ni - in_len + 1
                mj = zrodla(nrz)%polozenia(zrodla(nrz)%N,2)
                if(nj==1) nj = 2
                if(mj==NY)mj = NY-1
                GFLAGS(mi:ni,nj:mj) = B_NORMAL
            endif
        enddo

        ! -------------------------------------------------------------
		! wygladzanie krawedzi
        ! -------------------------------------------------------------
		do soft_iter = 1 , smooth_iter
        do i =  1+smooth_rad*2 , nx-smooth_rad*2-1
        do j =  1+smooth_rad , ny-smooth_rad-1
         norm_iter  = 0
		 empty_iter = 0
			do iii =  -smooth_rad , smooth_rad , 1
			do jjj =  -smooth_rad , smooth_rad , 1
				if(GFLAGS(i+iii,j+jjj) == B_EMPTY )  empty_iter = empty_iter + 1
				if(GFLAGS(i+iii,j+jjj) == B_NORMAL ) norm_iter = norm_iter   + 1
			enddo
			enddo
			if(empty_iter < norm_iter .and. GFLAGS(i,j) == B_EMPTY ) GFLAGS(i,j) = B_NORMAL
        enddo
        enddo
		enddo ! end of iter

		! tworzenie odpowiednich warunkow brzegowych (tak zeby nie bylo poloczenia
		! B_normal z B_Empty)
        do i = 2 , nx-1
        do j = 2 , ny-1
        		if (GFLAGS(i,j) == B_NORMAL) then
        			iter = (GFLAGS(i-1,j)-B_EMPTY) * (GFLAGS(i+1,j)-B_EMPTY) * (GFLAGS(i,j-1)-B_EMPTY) * (GFLAGS(i,j+1)-B_EMPTY)
        			if ( iter == 0 ) GFLAGS(i,j) = B_DIRICHLET
                endif
        enddo
        enddo

        ! usuwanie zbednych flag dirichleta
        do i = 1 , nx
        do j = 1 , ny
        		if (GFLAGS(i,j) == B_DIRICHLET) then
                    iter = 1
                    if(i>1 ) iter = iter*GFLAGS(i-1,j)*(GFLAGS(i-1,j)-B_WEJSCIE)
                    if(j>1 ) iter = iter*GFLAGS(i,j-1)*(GFLAGS(i,j-1)-B_WEJSCIE)
                    if(i<NX) iter = iter*GFLAGS(i+1,j)*(GFLAGS(i+1,j)-B_WEJSCIE)
                    if(j<NY) iter = iter*GFLAGS(i,j+1)*(GFLAGS(i,j+1)-B_WEJSCIE)
        			if ( iter /= 0 ) GFLAGS(i,j) = B_EMPTY
                endif
        enddo
        enddo

        ! nadawanie indeksow
        GINDEX = 0
        iter   = 1
        do i = 1 , nx
        do j = 1 , ny
           if( GFLAGS(i,j) == B_WEJSCIE .or. GFLAGS(i,j) == B_NORMAL .or. GFLAGS(i,j) == B_DIRICHLET  ) then
           GINDEX(i,j) = iter
           iter = iter + 1
           endif
        enddo
        enddo

    end subroutine system_inicjalizacja_ukladu


    ! --------------------------------------------------------------------
    ! Rozwiazywanie problemu dla zrodla o podanym numerze - nrz
    ! --------------------------------------------------------------------
    subroutine system_rozwiaz_problem(nrz)
        integer,intent(in)  :: nrz
        integer,allocatable :: HBROWS(:)
        integer :: i,j,i_comp,ni,nj,mod_in

        Ef  = atomic_Ef/1000.0/Rd
        BZ  = BtoDonorB(atomic_Bz)
        print*,"! ----------------------------------------------- !"
        print*,"! Rozpoczecie obliczen dla zrodla=", nrz
        print*,"!   Ef=", Ef*Rd*1000    , "[meV]"
        print*,"!   Bz=", DonorBtoB(Bz) , "[T]"
        print*,"! ----------------------------------------------- !"

        call oblicz_rozmiar_macierzy()

        allocate(CMATA (MATASIZE    ))
        allocate(IDXA  (MATASIZE,2  ))
        allocate(VPHI  (TRANS_MAXN  ))
        allocate(HBROWS(TRANS_MAXN+1))

        call system_generate_Wkl() ! wypelnianie macierzy Wkl

        mod_in = 1
        zrodla(nrz)%ck(:) = 0
        zrodla(nrz)%ck(mod_in) = 1
        call zrodla(nrz)%zrodlo_oblicz_Fj(dx)
        VPHI = 0
        do i = 1 , zrodla(nrz)%N
            ni = zrodla(nrz)%polozenia(i,1)
            nj = zrodla(nrz)%polozenia(i,2)
            VPHI(GINDEX(ni,nj)) = zrodla(nrz)%Fj(i)
        enddo
        call wypelnij_macierz()
        call convert_to_HB(MATASIZE,IDXA,HBROWS)
        call solve_system(TRANS_MAXN,MATASIZE,IDXA(:,2),HBROWS,CMATA(:),VPHI)

        call oblicz_TR(nrz,mod_in)

        PHI  = 0
        do i = 1 , Nx
        do j = 1 , Ny
            do i_comp = 1 , no_components
                if(GINDEX(i,j) > 0) PHI(i,j,i_comp) = VPHI(GINDEX(i,j)+(i_comp-1)*TRANS_COMPOMENT_MAXN)
            enddo
        enddo
        enddo


        deallocate(CMATA)
        deallocate(IDXA)
        deallocate(VPHI)
        deallocate(HBROWS)


    end subroutine system_rozwiaz_problem

    ! ---------------------------------------------------------------------
    !   Wyliczanie rozmiaru macierzy: zakladamy, ze na :
    ! 1. warunek Dirichleta potrzeba jedna wartosc
    ! 2. normalny punkt wewnatrz ukladu potrzeba 5 wartosci (4 siasiadow + 1 wezel centralny)
    ! 3. dla wejsc potrzeba N-2 + 1 wezlow: gdzie N-2 oznacza ze bierzemy
    !    wezly wejsciowe od 2 .. N -1, a +1 wynika z faktu, ze mamy zawsze
    !    polaczenie z wezlem (i+1,j), lub (i-1,j)
    ! ---------------------------------------------------------------------

    subroutine oblicz_rozmiar_macierzy()
        integer :: i , j
        MATASIZE   = 0
        TRANS_MAXN = 0
        do i = 1 , NX
        do j = 1 , Ny
            select case(GFLAGS(i,j))
            case(B_DIRICHLET)
                    MATASIZE   = MATASIZE   + 1
                    TRANS_MAXN = TRANS_MAXN + 1
            case(B_WEJSCIE)
                    MATASIZE   = MATASIZE   + zrodla(ZFLAGS(i,j))%N + 1 - 2 ! +1 na relacje prawo/lewo -2 na obciecie indeksow (1 i N)
                    TRANS_MAXN = TRANS_MAXN + 1
            case(B_NORMAL)
                    MATASIZE   = MATASIZE   + 5
                    TRANS_MAXN = TRANS_MAXN + 1
            end select
        enddo
        enddo
        ! przeliczamy rozmiary macierzy na nowe z uwzglednieniem tego, ze
        ! mamy kilka nowych skladowych f.f. : rozmiar macierzy rosnie no_components-razy plus
        ! elementy diagonalne (W_kl)
        MATASIZE   = MATASIZE * no_components  + (no_components**2 - no_components)*TRANS_MAXN
        ! rozmiar f.f zwieksza sie no_components razy
        TRANS_COMPOMENT_MAXN = TRANS_MAXN
        TRANS_MAXN = TRANS_MAXN * no_components

        print*,"Rozmiar macierzy - H:",MATASIZE
        print*,"Rozmiar problemu - Y:",TRANS_MAXN
        print*,"Rozmiar komponentu-Y:",TRANS_COMPOMENT_MAXN
    end subroutine oblicz_rozmiar_macierzy


    ! ---------------------------------------------------------------------
    !   Funkcja wypelnia macierz hamiltonianu w sposob rzadki, zapisujac
    ! kolejne niezerowe elementy H do dwoch wektorow:
    ! cmatA(i) - i-ta wartosc niezerowa w macierzy
    ! idxA(i,1) - w jakim wierszu macierzy sie znajduje i-ty nie zerowy element
    ! idxA(i,2) - oraz w jakiej kolumnie
    ! W zaleznosci od flagi wezla elementy macierzy H sa wyznaczane w
    ! odpowiedni sposob.
    ! ---------------------------------------------------------------------

    subroutine wypelnij_macierz()

        ! zmienne pomocniczne
        integer          :: i,j,itmp,ni,nj,pnj,nn,ln,nzrd
        integer          :: component_iter , component_index_offset , col_component_iter
        complex*16       :: post
        itmp  = 1

        component_index_offset = TRANS_COMPOMENT_MAXN ! przesuniecie w indeksach
        do component_iter = 1 , no_components

        do i = 1, nx
        do j = 1, ny
            if(GINDEX(i,j) > 0) then
                if( GFLAGS(i,j) == B_DIRICHLET ) then

                    cmatA(itmp)  = CMPLX(1.0)
                    idxA(itmp,1) = GINDEX(i, j) + (component_index_offset)*(component_iter-1)
                    idxA(itmp,2) = GINDEX(i, j) + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    ! dodajemy elementy diagonalne w macierzy ale tylko dla roznych skladowych W_kl
                    do col_component_iter =  1 , no_components
                        if( component_iter /= col_component_iter ) then ! sprawdzamy czy nie jestesmy na diagonali
                            cmatA(itmp)  = W_kl(component_iter,col_component_iter,GINDEX(i, j))
                            idxA(itmp,1) = GINDEX(i, j) + (component_index_offset)*(component_iter    -1)
                            idxA(itmp,2) = GINDEX(i, j) + (component_index_offset)*(col_component_iter-1)
                            itmp = itmp + 1
                        endif
                    enddo

                elseif( GFLAGS(i,j) == B_NORMAL) then
                    cmatA(itmp)   = CMPLX( 2.0/DX/DX + UTOTAL(i,j) - Ef )
                    idxA (itmp,1) = GINDEX(i,j) + (component_index_offset)*(component_iter-1)
                    idxA (itmp,2) = GINDEX(i,j) + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    cmatA(itmp)   = CMPLX(-0.5/DX/DX*EXP(+II*DX*DX*(j)*BZ)  )
                    idxA (itmp,1) = GINDEX(i  ,j) + (component_index_offset)*(component_iter-1)
                    idxA (itmp,2) = GINDEX(i-1,j) + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    cmatA(itmp)   = CMPLX(-0.5/DX/DX*EXP(-II*DX*DX*(j)*BZ)  )
                    idxA (itmp,1) = GINDEX(i  ,j) + (component_index_offset)*(component_iter-1)
                    idxA (itmp,2) = GINDEX(i+1,j) + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    cmatA(itmp) = CMPLX(-0.5/DX/DX)
                    idxA(itmp,1) = GINDEX(i,j)    + (component_index_offset)*(component_iter-1)
                    idxA(itmp,2) = GINDEX(i,j+1)  + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    cmatA(itmp) = CMPLX(-0.5/DX/DX)
                    idxA(itmp,1) = GINDEX(i,j)    + (component_index_offset)*(component_iter-1)
                    idxA(itmp,2) = GINDEX(i,j-1)  + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    ! dodajemy elementy diagonalne w macierzy ale tylko dla roznych skladowych W_kl
                    do col_component_iter =  1 , no_components
                        if( component_iter /= col_component_iter ) then ! sprawdzamy czy nie jestesmy na diagonali
                            cmatA(itmp)  = W_kl(component_iter,col_component_iter,GINDEX(i, j))
                            idxA(itmp,1) = GINDEX(i, j) + (component_index_offset)*(component_iter    -1)
                            idxA(itmp,2) = GINDEX(i, j) + (component_index_offset)*(col_component_iter-1)
                            itmp = itmp + 1
                        endif
                    enddo
                ! ----------------------------------------------------------------------
                ! Obsluga wejsc
                ! ----------------------------------------------------------------------
                else if( GFLAGS(i,j) == B_WEJSCIE) then
                    nzrd = ZFLAGS(i,j)
                    ! ------------------------------------------------------------------
                    ! Zrodla prawe:
                    !
                    !                      ------------------------>
                    !
                    ! ------------------------------------------------------------------
                    if( zrodla(nzrd)%bKierunek == .true. ) then

                    ni   = i
                    nj   = j ! globalne polozenia na siatce
                    ln   = GINDEX(i,j) - GINDEX(i,zrodla(nzrd)%polozenia(1,2)) + 1 ! lokalny indeks

                    cmatA(itmp)   = CMPLX(-0.5/DX/DX*2*cos(DX*DX*(nj)*BZ))
                    idxA (itmp,1) = GINDEX(i  ,j) + (component_index_offset)*(component_iter-1)
                    idxA (itmp,2) = GINDEX(i+1,j) + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    post = 0.5/DX/DX*(EXP(II*DX*DX*(nj)*BZ))

                    do nn = 2 , zrodla(nzrd)%N - 1
                    pnj   = zrodla(nzrd)%polozenia(nn,2)
                    if( ln == nn  ) then

                        cmatA(itmp)   = CMPLX( 2.0/DX/DX + UTOTAL(i,j) - Ef ) &
                           & + post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA (itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA (itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1

                    else if( nn == ln+1  ) then
                        cmatA(itmp) = CMPLX(-0.5/DX/DX) &
                        &   + post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA(itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA(itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1
                    else if( nn == ln-1  ) then

                        cmatA(itmp) = CMPLX(-0.5/DX/DX) &
                        &   + post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA(itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA(itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1
                    else
                        cmatA(itmp)  = post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA(itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA(itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1
                    endif
                    enddo

                    ! dodajemy elementy diagonalne w macierzy ale tylko dla roznych skladowych W_kl
                    do col_component_iter =  1 , no_components
                        if( component_iter /= col_component_iter ) then ! sprawdzamy czy nie jestesmy na diagonali
                            cmatA(itmp)  = W_kl(component_iter,col_component_iter,GINDEX(i, j))
                            idxA(itmp,1) = GINDEX(i, j) + (component_index_offset)*(component_iter    -1)
                            idxA(itmp,2) = GINDEX(i, j) + (component_index_offset)*(col_component_iter-1)
                            itmp = itmp + 1
                        endif
                    enddo

                    else
                    ! ------------------------------------------------------------------
                    ! Zrodla lewe:
                    !
                    !                      <------------------------
                    !
                    ! ------------------------------------------------------------------
                    ni   = i
                    nj   = j ! globalne polozenia na siatce
                    ln   = GINDEX(i,j) - GINDEX(i,zrodla(nzrd)%polozenia(1,2)) + 1 ! lokalny indeks

                    cmatA(itmp)   = CMPLX(-0.5/DX/DX*2*cos(DX*DX*(nj)*BZ))
                    idxA (itmp,1) = GINDEX(i  ,j) + (component_index_offset)*(component_iter-1)
                    idxA (itmp,2) = GINDEX(i-1,j) + (component_index_offset)*(component_iter-1)
                    itmp = itmp + 1

                    post =-0.5/DX/DX*(EXP(-II*DX*DX*(nj)*BZ))

                    do nn = 2 , zrodla(nzrd)%N - 1
                    pnj   = zrodla(nzrd)%polozenia(nn,2)
                    if( ln == nn  ) then

                        cmatA(itmp)   = CMPLX( 2.0/DX/DX + UTOTAL(i,j) - Ef ) &
                           & + post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA (itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA (itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1

                    else if( nn == ln+1  ) then
                        cmatA(itmp) = CMPLX(-0.5/DX/DX) &
                        &   + post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA(itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA(itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1
                    else if( nn == ln-1  ) then

                        cmatA(itmp) = CMPLX(-0.5/DX/DX) &
                        &   + post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA(itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA(itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1
                    else
                        cmatA(itmp)  = post*zrodla(nzrd)%zrodlo_alfa_v_i(dx,ln,nn)
                        idxA(itmp,1) = GINDEX(ni,nj)  + (component_index_offset)*(component_iter-1)
                        idxA(itmp,2) = GINDEX(ni,pnj) + (component_index_offset)*(component_iter-1)
                        itmp = itmp + 1
                    endif
                    enddo

                    ! dodajemy elementy diagonalne w macierzy ale tylko dla roznych skladowych W_kl
                    do col_component_iter =  1 , no_components
                        if( component_iter /= col_component_iter ) then ! sprawdzamy czy nie jestesmy na diagonali
                            cmatA(itmp)  = W_kl(component_iter,col_component_iter,GINDEX(i, j))
                            idxA(itmp,1) = GINDEX(i, j) + (component_index_offset)*(component_iter    -1)
                            idxA(itmp,2) = GINDEX(i, j) + (component_index_offset)*(col_component_iter-1)
                            itmp = itmp + 1
                        endif
                    enddo

                    endif ! end of if PRAWE/ LEWE
                endif ! end of if WEJSCIE
            endif ! end if GINDEX > 0
        enddo ! end of j
        enddo ! end of i

        end do ! end of component_iter - petla po skladowych f.f.

        if ( itmp -1 - MATASIZE /= 0 ) then
            print *, "Blad rozmiar maicierzy niezgodny z przewidywanym:" , itmp - 1
            stop
        else
            print *, "Macierz skonstruowana poprawnie..."
        endif


    end subroutine wypelnij_macierz

    ! -------------------------------------------------------------------
    ! Procedura oblicza prawdopodobienstwo przejscia T oraz odbicia R
    ! dla podanego zrodla wejsciowego nrz i danego modu wchodzacego mod_in
    ! -------------------------------------------------------------------
    subroutine oblicz_TR(nrz,mod_in)
        integer,intent(in) :: nrz,mod_in
        integer :: i,i_comp
        doubleprecision :: Jin


        if( TRANS_DEBUG == .true. )  print*,"obliczenia dla pierwszej skladowej:"
        ! obliczamy apliduty prawdopodobienstwa
        do i = 1 , no_zrodel
            call zrodla(i)%zrodlo_oblicz_dk(VPHI,GINDEX,dx)
            !call zrodla(i)%zrodlo_wypisz_ckdk()
        enddo
        ! na podstawie amplitud prawdopodobienstwa obliczamy strumienie
        do i = 1 , no_zrodel
            call zrodla(i)%zrodlo_oblicz_JinJout(dx)
        enddo

        ! na podstawie strumieni obliczamy transmisje oraz prawdopodobienstwo odbicia
        Jin = zrodla(nrz)%Jin(mod_in)
        do i = 1 , no_zrodel
            zrodla(i)%Jin(:)  = zrodla(i)%Jin(:)   / Jin
            zrodla(i)%Jout(:) = zrodla(i)%Jout(:)  / Jin
            call zrodla(i)%zrodlo_wypisz_JinJout()
        enddo

        TRANS_R = sum( abs(zrodla(nrz)%Jout(:)) )
        TRANS_T = 0
        do i = 1 ,no_zrodel
            if(i /= nrz) TRANS_T = TRANS_T + sum( abs(zrodla(i)%Jout(:)) )
        enddo
        TRANS_R_COMPONENTS = 0
        TRANS_T_COMPONENTS = 0
        TRANS_R_COMPONENTS(1) = TRANS_R
        TRANS_T_COMPONENTS(1) = TRANS_T




        zrodla(nrz)%ck(:) = 0
        if( TRANS_DEBUG == .true. ) then
            print*,"T[1]= ", TRANS_T
            print*,"R[1]= ", TRANS_R
            print*,"W[1]= ", TRANS_R + TRANS_T
            print*,"Obliczenia dla drugiej skladowej:"
        endif
        do i_comp = 2 , no_components

            ! obliczamy apliduty prawdopodobienstwa
            do i = 1 , no_zrodel
                call zrodla(i)%zrodlo_oblicz_dk(VPHI,GINDEX + (i_comp-1)*TRANS_COMPOMENT_MAXN,dx)
                !call zrodla(i)%zrodlo_wypisz_ckdk()
            enddo
            ! na podstawie amplitud prawdopodobienstwa obliczamy strumienie
            do i = 1 , no_zrodel
                call zrodla(i)%zrodlo_oblicz_JinJout(dx)
            enddo

            ! na podstawie strumieni obliczamy transmisje oraz prawdopodobienstwo odbicia
            do i = 1 , no_zrodel
                zrodla(i)%Jin(:)  = zrodla(i)%Jin(:)   / Jin
                zrodla(i)%Jout(:) = zrodla(i)%Jout(:)  / Jin
                call zrodla(i)%zrodlo_wypisz_JinJout()
            enddo

            TRANS_R = 0
            TRANS_T = 0
            do i = 1 ,no_zrodel
                if(      zrodla(i)%bKierunek) TRANS_T = TRANS_T + sum( abs(zrodla(i)%Jout(:)) )
                if(.not. zrodla(i)%bKierunek) TRANS_R = TRANS_R + sum( abs(zrodla(i)%Jout(:)) )
            enddo
            TRANS_R_COMPONENTS(i_comp) = TRANS_R
            TRANS_T_COMPONENTS(i_comp) = TRANS_T

            if( TRANS_DEBUG == .true. ) then
                print*,"T[",i_comp,"]= ", TRANS_T
                print*,"R[",i_comp,"]= ", TRANS_R
                print*,"W[",i_comp,"]= ", TRANS_R + TRANS_T
            endif

        enddo ! end of i_comp loop

        TRANS_T = sum(TRANS_T_COMPONENTS)
        TRANS_R = sum(TRANS_R_COMPONENTS)
        print*,"T = ", TRANS_T
        print*,"R = ", TRANS_R
        print*,"W = ", TRANS_R + TRANS_T

    end subroutine oblicz_TR
    ! -------------------------------------------------------------------
    ! Zwraca wartosc rozkladu gaussa w punkcie (x,y) dla gaussa o srodku
    ! w (xpos,ypos) oraz sigmie = sigma i amplidudzie = amplitude.
    ! x,y,xpos,ypos - podajemy w indeksach siatki (i,j)
    ! -------------------------------------------------------------------
    doubleprecision function system_gauss(x,y,xpos,ypos,sigma,amplitude)
        integer :: x,y,xpos,ypos
        doubleprecision :: sigma,amplitude
            system_gauss = amplitude*exp(-sigma*(( x - xpos )**2+( y - ypos )**2)*DX*DX)
    end function system_gauss

    ! --------------------------------------------------------------------
    !
    ! --------------------------------------------------------------------
    subroutine system_generate_Wkl()
            integer :: i,j,ix,iy,k , no_nbr
            integer :: xpos , ypos , norand
            doubleprecision :: hval , gamp
            complex*16 :: ave_wkl
            W_kl = 0

            do i = 1 , no_components
            do j = 1 , no_components
!                do norand = 1 , 120
!                xpos =  rand() * (Nx-40) + 20
!                ypos =  rand() * Ny + 1
!                hval =  40.0
!                gamp = -1/dx/dx*(rand()-0.5)*0.3
!
!                do ix = 1 , NX
!                do iy = 1 , NY
!                    if( ix > 10 .and. ix < nx-10 ) then
!                        if (GINDEX(ix,iy) > 0 ) W_kl(i,j,GINDEX(ix,iy)) = W_kl(i,j,GINDEX(ix,iy)) + system_gauss(ix,iy,xpos,ypos,hval,gamp)
!                        !if (GINDEX(ix,iy) > 0 ) W_kl(i,j,GINDEX(ix,iy)) = 1.0*1/dx/dx*(rand()-0.5)
!                    endif
!                enddo
!                enddo
!                enddo

                do ix = 1 , NX
                do iy = 1 , NY
                    if( ix > nx/10 .and. ix < nx-nx/10 ) then
                        !if (GINDEX(ix,iy) > 0 ) W_kl(i,j,GINDEX(ix,iy)) = W_kl(i,j,GINDEX(ix,iy)) + system_gauss(ix,iy,xpos,ypos,hval,gamp)
                        if (GINDEX(ix,iy) > 0 ) W_kl(i,j,GINDEX(ix,iy)) = decoh_amp*1/dx/dx*(rand()-0.5)
                    endif
                enddo
                enddo
            enddo
            enddo
            do k = 1, 0 ! sredniowanie
            do i = 1 , no_components
            do j = 1 , no_components
                do ix = 2 , NX-1
                do iy = 1 , NY
                    if( ix > 10 .and. ix < nx-10 ) then
                        ave_wkl = 0
                        no_nbr = 0
                        if (GINDEX(ix+1,iy) > 0 ) then
                            ave_wkl = ave_wkl + W_kl(i,j,GINDEX(ix+1,iy))
                            no_nbr = no_nbr + 1
                            endif
                        if (GINDEX(ix-1,iy) > 0 ) then
                            ave_wkl = ave_wkl + W_kl(i,j,GINDEX(ix-1,iy))
                            no_nbr = no_nbr + 1
                            endif
                        if(iy<NY) then
                            if (GINDEX(ix,iy+1) > 0 ) then
                            ave_wkl = ave_wkl + W_kl(i,j,GINDEX(ix,iy+1))
                            no_nbr = no_nbr + 1
                            endif
                            endif
                        if(iy> 1) then
                            if (GINDEX(ix,iy-1) > 0 ) then
                            ave_wkl = ave_wkl + W_kl(i,j,GINDEX(ix,iy-1))
                            no_nbr = no_nbr + 1
                            endif
                            endif
                        if (GINDEX(ix,iy)   > 0 ) W_kl(i,j,GINDEX(ix,iy))  = ave_wkl/no_nbr
                    endif
                enddo
                enddo
            enddo
            enddo
            enddo ! end of k

            ! symetryzacja do macierzy hermitowskiej
            do i = 1 , no_components
            do j = i , no_components
                do k = 1 , TRANS_COMPOMENT_MAXN
                    W_kl(i,j,k) = conjg(W_kl(j,i,k))
                enddo
            enddo
            enddo
    end subroutine system_generate_Wkl


    ! ------------------------------------------------------------ -------
    ! Funkcja zapisuje do pliku nazwa dane wskazane przez typ (patrz enum)
    ! ------------------------------------------------------------ -------
    subroutine system_zapisz_do_pliku(nazwa,typ)
    character(*) , intent(in) :: nazwa
    integer                   :: typ

    double precision :: fval
    double precision :: comp_val(no_components)
    complex*16       :: wkl_val(no_components,no_components)
    integer          :: i,j,i_comp

    open(unit=86554,file=nazwa)
    select case(typ)
    ! ---------------------------------------
    case(ZAPISZ_KONTUR)
        do i = 1 , NX
        do j = 1 , NY
                fval = 0;
                if(GINDEX(i,j) > 0) fval = 1
                write(86554,"(3e20.6)"),i*DX*Lr2L,j*DX*Lr2L,fval
        enddo
            write(86554,*),""
        enddo
    ! ---------------------------------------
    case(ZAPISZ_FLAGI)
        do i = 1 , NX
        do j = 1 , NY
                fval = GFLAGS(i,j)
                write(86554,"(3e20.6)"),i*DX*Lr2L,j*DX*Lr2L,fval
        enddo
            write(86554,*),""
        enddo
    ! ---------------------------------------
    case(ZAPISZ_POTENCJAL)
        do i = 1 , NX
        do j = 1 , NY
                fval = UTOTAL(i,j)
                write(86554,"(3e20.6)"),i*DX*Lr2L,j*DX*Lr2L,fval*Rd*1000
        enddo
            write(86554,*),""
        enddo
    ! ---------------------------------------
    case(ZAPISZ_INDEKSY)
        do i = 1 , NX
        do j = 1 , NY
                fval = GINDEX(i,j)
                write(86554,"(3e20.6)"),i*DX*Lr2L,j*DX*Lr2L,fval
        enddo
            write(86554,*),""
        enddo
    ! ---------------------------------------
    case(ZAPISZ_PHI)
        do i = 1 , NX
        do j = 1 , NY
                do i_comp = 1 , no_components
                comp_val(i_comp) = abs(PHI(i,j,i_comp))**2
                enddo
                write(86554,"(100e20.6)"),i*DX*Lr2L,j*DX*Lr2L,comp_val
        enddo
            write(86554,*),""
        enddo
    case(ZAPISZ_WKL)
        do i = 1 , NX
        do j = 1 , NY
                wkl_val = 0
                if( GINDEX(i,j) > 0 ) wkl_val = W_kl(:,:,GINDEX(i,j))
                write(86554,"(100f20.6)"),i*DX*Lr2L,j*DX*Lr2L,dble(wkl_val)*Rd*1000
        enddo
            write(86554,*),""
        enddo
    ! ---------------------------------------
    case default
            print*,"System: Zapisz do pliku - podano zly argument."
    end select

    close(86554)
    endsubroutine system_zapisz_do_pliku

    ! ==========================================================================
    !
    !
    !                          SUPER LU
    !
    !
    !
    ! ==========================================================================

    subroutine convert_to_HB(no_vals,rows_cols,out_rows)
          integer,intent(in)                  :: no_vals
          integer,intent(in),dimension(:,:)  :: rows_cols
          integer,intent(inout),dimension(:) :: out_rows
          integer :: iterator, irow
          integer :: i, n

          n        = no_vals
          iterator = 0
          irow     = 0
          do i = 1 , n
              if( rows_cols(i,1) /= irow ) then
                iterator = iterator + 1
                out_rows(iterator) = i
                irow = rows_cols(i,1)
              endif
          enddo
          out_rows(iterator+1) = n + 1
      end subroutine convert_to_HB


    subroutine solve_system(no_rows,no_vals,colptr,rowind,values,b)
        integer,intent(in)                 :: no_rows
        integer,intent(in)                 :: no_vals
        integer,intent(in),dimension(:)    :: colptr,rowind
        complex*16,intent(in),dimension(:) :: values
        complex*16,intent(inout),dimension(:) :: b

        integer n, nnz, nrhs, ldb, info, iopt
        integer*8 factors
!
!      call zhbcode1(n, n, nnz, values, rowind, colptr)
!

        n    = no_rows
        nnz  = no_vals
        ldb  = n
        nrhs = 1


! First, factorize the matrix. The factors are stored in *factors* handle.
      iopt = 1
      call c_fortran_zgssv( iopt, n, nnz, nrhs, values, colptr , rowind , b, ldb,factors, info )
!
      if (info .eq. 0) then
         write (*,*) 'Factorization succeeded'
      else
         write(*,*) 'INFO from factorization = ', info
      endif
!
! Second, solve the system using the existing factors.
      iopt = 2
      call c_fortran_zgssv( iopt, n, nnz, nrhs, values, colptr,rowind ,  b, ldb,factors, info )
!
      if (info .eq. 0) then
         write (*,*) 'Solve succeeded'
!         write (*,*) (b(i), i=1, n)
      else
         write(*,*) 'INFO from triangular solve = ', info
      endif

! Last, free the storage allocated inside SuperLU
      iopt = 3
      call c_fortran_zgssv( iopt, n, nnz, nrhs, values, colptr,rowind, b, ldb,factors, info )

      endsubroutine solve_system

end module
